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Abstract 

This  report  presents  a  computational  technique  for  the  rapid,  efficient  calculation  of 
fields  from  transient  acoustic  sources  in  linear,  isotropic  media.  The  source  velocity  is 
separable  in  space  and  time.  The  method  uses  a  spatial  impulse  response  method  based  on 
linear  systems  concepts  to  express  the  output  in  terms  of  the  Green's  function  of  propagation 
equation  and  the  boundary  conditions.  The  output  is  expressed  as  the  inverse  spatial 
transform  of  the  product  of  the  transform  of  the  spatial  excitation  and  a  time-varying 
spatial  filter  that  represents  propagation.  The  calculation  technique  has  been  implemented 
in  MATLAB  and  sample  cases  are  presented  for  the  circular  and  square  piston,  as  well  as  a 
Gaussian-  and  Bessel-weighted  spatial  excitation.  Code  for  the  MATLAB  implementation 
is  provided. 
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Chapter  1 

Introduction 

1.1      Introduction 

The  propagation  characteristics  of  continuously  radiated  monochromatic  ultrasonic  sources 
are  well  solved  through  application  of  the  angular  spectrum  technique  [l]  or  Fresnel  in- 
tegrals. More  frequently,  however,  acoustic  imaging,  tissue  characterization,  and  physical 
acoustics  applications  tend  to  use  pulsed  sound.  The  propagation  of  pulsed  ultrasound  with 
arbitrary  temporal  and  spatial  components  is  not  understood  to  the  same  degree.  We  want 
to  develop  a  reliable,  easy  method  of  diffraction  prediction.  This  report  describes  an  ap- 
proach based  on  linear  systems  theory  and  the  Fourier  transform.  The  goal  was  to  achieve  a 
readily  usable  method  of  predicting  pulsed  wave  diffraction  in  a  time-efficient  and  accurate 
manner  in  order  to  examine  the  wave  diffraction.  Earlier  efforts  had  used  FORTRAN  code 
to  implement  the  propagation  model  on  both  a  mainframe  computer  [2]-  [6]  and  a  personal 
computer  [7].  Merrill  had  attempted  to  use  MATLAB  (a  commercial  matrix  manipulation 
software  package)  to  perform  the  calculations  on  a  PC  but  was  thwarted  by  memory  re- 
strictions of  the  earlier  versions  of  this  software.  More  recent  versions  allow  larger  matrices 
and  arrays  to  be  manipulated  (subject  only  to  the  computer  memory  available)  Thesis  work 
performed  by  Upton  [8]  and  Reid  [9]  applied  this  updated  software  to  the  problems  of  op- 
tical propagation  and  acoustic  propagation,  respectively.  Much  of  the  foundation  for  the 
work  reported  here  was  done  by  these  officer-students. 

The  basic  method  of  the  spatial  impulse  response  was  introduced  by  Stepanishen  [10]- 
[13],  reviewed  by  Harris  [14],  and  modified  by  Guyomar  and  Powers  [5].  The  approach  of 
Guyomar  and  Powers  differed  by  the  use  of  linear  systems  theory.  Linear  systems  theory 
revealed  the  importance  of  the  total  impulse  response  and  its  equivalence  to  the  Green's 
function.  Furthermore,  the  spatial  impulse  response  functions  are  found  in  the  spatial 
transform  (Fourier)  domain.  Working  in  the  transform  domain  allowed  propagation  of  the 
wave  to  be  viewed  as  a  time-varying  spatial  filter  applied  to  the  spatial  spectrum  of  the 
input  excitation.  The  advantage  of  this  method  is  that  it  provides  the  diffracted  field  from 
an  insonifying  wave  with  arbitrary  temporal  and  spatial  dependence  in  a  computationally 


efficient  form.  By  use  of  the  Fourier  transform,  an  efficient  computer  implementation  of 
this  technique  using  FFT  routines  is  possible. 

The  desired  benefit  of  a  fast,  time-efficient  computer  implementation  to  calculate  the 
acoustic  potential  or  pressure  is  to  aid  in  ultrasonic  transducer  design  for  medical,  acoustic 
imaging,  and  mine  warfare  applications.  With  a  knowledge  of  wave  diffraction  phenomenon 
a  diffracted  wave  reflected  from  an  unknown  object  can  be  used  to  provide  information 
about  the  object.  This  type  of  system  must  be  portable  as  well  as  time-efficient,  which  is 
very  achievable  given  the  trends  in  computer  technology.  Computers  have  become  faster  and 
have  increased  memory  capacity  while  their  size  has  decreased.  Other  benefits  are  derived 
from  the  use  of  the  matrix  manipulation  program,  MATLAB,  and  the  ability  to  expand 
this  implementation  to  cases  involving  lossy  media.  Because  MATLAB  is  readily  available 
on  the  commercial  market,  it  requires  no  special  equipment  for  computer  implementation. 

1.2      Report  Overview 

Chapter  II  consists  of  the  problem  description,  including  the  source-to-observation  plane 
geometry,  a  discussion  on  the  linear  systems  approach,  the  mathematical  development,  and 
an  overview  of  MATLAB  and  AXUM.  Chapter  III  consists  of  a  discussion  of  the  two  pro- 
gram modules,  the  acoustic  filter  module,  AC-FIL.M,  and  the  acoustic  propagation  module, 
AC-PROP.M.  Chapter  IV  starts  with  the  set  of  defining  parameters.  These  parameters  are 
then  used  for  an  explanation  of  the  program's  verification  and  an  investigation  of  other 
input  excitation  functions.  Chapter  V  contains  some  examples  of  the  numerical  output  for 
various  source  excitation  conditions. 

Following  a  summary  in  Chapter  VI,  Appendices  A  and  C  give  detailed  explanations 
of  the  MATLAB  routines  that  were  used  to  produce  the  numerical  results,  AC-FIL.M  and 
ACPROP.M.  Appendix  B  gives  examples  of  the  propagation  filters  generated  by  the  code 
in  Appendix  A.  The  source  code  of  the  various  geometric  excitation  functions  is  given  in 
Appendix  D.  Representative  output  presentations  using  a  scientific  visualization  program, 
called  Spyglass  Dicer,  are  found  in  Appendix  E. 


Chapter  2 

Problem  Description 


Before  assembling  a  computer  implementation,  we  must  first  understand  the  problem.  An 
explanation  of  the  problem  follows,  beginning  with  the  description  of  the  geometry  in  the 
first  section.  Section  two  continues  the  explanation  into  the  linear  systems  approach.  The 
third  section  proceeds  through  the  mathematical  development  of  the  problem  and  ties  in 
the  Fourier  transform.  The  theory  presented  in  the  first  three  sections  was  derived  from 
the  works  of  Guyomar  and  Powers  [3,  4,  5,  6].  The  final  section  gives  an  overview  of  the 
software  tools  used  for  generating  the  excitation  functions,  the  propagation  fields,  and  the 
output  graphics. 

2.1  Geometry 

The  problem  geometry  is  shown  in  Fig.  2.1.  The  acoustic  velocity  potential  <p(x,y,z)  is 
to  be  calculated  at  an  arbitrary  point  in  the  positive- z  half-space  given  the  z-directed 
velocity  in  a  portion  of  the  source  plane.  The  source's  z-directed  velocity  is  spatially  and 
temporally  arbitrary  and  is  rigidly  baffled  (i.e.,  equal  to  zero)  in  the  region  outside  the 
source.  Furthermore,  it  is  assumed  that  the  spatial  and  temporal  components  of  the  z- 
velocity  are  separable,  having  the  form  at  the  input  plane 

v2{xo, i/o,0,0    =     T(t)s(x0,y0)  •  (2.1) 

A  linear,  homogeneous,  and  lossless  medium  is  assumed  to  be  present  between  the  source 
and  observation  point.  (The  extension  to  lossy  propagation  can  be  found  in  Ref.  [15].) 

2.2  Linear  Systems  Approach 

Linear  systems  solutions  are  applied  to  systems  that  are  linear  and  time-invariant.  A 
linear  systems  solution  approach  to  this  problem  is  possible  because  propagation  in  a  linear 
homogeneous  media  is  a  linear,  space-invariant  process  [3].  In  linear  systems  the  impulse 
response  is  the  response  of  the  system  to  an  impulsive  input.   The  total  impulse  response 


Figure  2.1:  Source-to-receiver  geometry. 


h(x,y,z,t)  of  a  system  is  produced  by  an  input  of  the  form  6(x,y,t)  =  6(x,y)6(t);  this  is 
shown  in  Fig.  2.2a.  Figure  2.2b  shows  the  spatial  impulse  response  p(x,y,z,t),  defined  as 
the  response  to  an  excitation  of  the  form  s(x,y)6(t).  Note  that  an  arbitrary  spatial  input 
has  been  substituted  for  the  impulsive  spatial  input.  Recall  from  linear  systems  theory  that 
the  solution  for  an  arbitrary  input  (spatial  or  temporal)  is  the  convolution  of  the  input  with 
the  system's  total  impulse  response;  therefore,  the  spatial  impulse  response  in  this  case  is 
given  by 


p(x,y,z,t)    =     s(x,y)mxmh(x,y,z,t), 


(2.2) 


where  *  indicates  convolution  with  respect  to  the  variable  shown. 

The  double  convolution  in  Eq.  2.2  is  not  easily  computer  implemented.  The  Fourier 
transform  furnishes  a  convenient  method  to  resolve  this  dilemma  by  using  the  property 
that  convolution  in  the  spatial  domain  becomes  multiplication  in  the  transform  domain, 
i.e., 


P{fxJy,Z,t)      =      s{fx,fy)h(fx,fy,Z,t), 


(2.3) 


where  the  tilde  notation  indicates  the  Fourier  transform  of  the  function  (in  this  case,  the 
two-dimensional  spatial  Fourier  transform).  This  relation  is  shown  in  Fig.  2.2c. 
The  general  solution  is  shown  in  Fig.  2.2d,  where 


<j>{x,y,z,t)     -     s(x,y)T(t)xmmth(x,y,z,t). 


(2.4) 
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Figure  2.2:  Block  diagram  of  (a)  the  total  impulse  response  h(x,y,t)),  (b)  the  spatial 
impulse  response  p(x,y,z,t)  in  the  space-time  domain,  (c)  the  spatial  impulse  response 
p(fx,fy,t)  in  the  spatial  frequency-time  domain,  and  (d)  the  general  solution  <p(x,y,z,t). 


Substituting  Eq.  2.2  into  Eq.  2.4  gives 

<t>(x,y,z,t)     =     T(t)'tp(x,y,z,t).  (2.5) 

It  follows  from  Fig.  2.2c  and  Eqs.  2.2  and  2.5  that  the  key  to  finding  the  general  solution 
is  finding  the  spatial  impulse  response  p(x,  y,z,t)  which  is,  in  turn,  dependent  on  the  total 
impulse  response  of  the  system  h(x,y,z,t).  This  total  impulse  response  of  the  system  is 
the  propagation  field  that  results  from  an  impulsive  source,  as  in  Fig  2.2a,  that  solves  the 
wave  equation  and  satisfies  the  boundary  conditions.  The  solution  to  the  wave  equation 
satisfying  the  boundary  conditions  is  commonly  known  as  the  Green's  function;  hence,  the 
total  impulse  response  is  simply  the  Green's  function.  Therefore,  once  the  Green's  function 
is  known,  the  total  impulse  response  is  also  known,  and  the  solution  becomes  a  triple 
convolution  between  an  excitation  source  which  is  spatially  and  temporally  arbitrary  (and 
assumed  to  be  known),  and  the  systems'  total  impulse  response  or  Green's  function.  The 
two  spatial  convolutions  are  easily  implemented  in  the  spatial  transform  domain;  the  time 
convolution  can  be  implemented  in  the  temporal  transform  domain,  if  desired.  We  have 
found  it  easier  to  perform  the  time  convolution  directly  in  the  time  domain. 

2.3      Mathematical  Development  for  Acoustic  Propagation 

The  wave  equation  for  lossless  media  is 

9  1  d2<t> 

The  general  solution  of  Fig.  2.2c  gives  the  result  in  terms  of  the  acoustic  potential  <p  which 
must  be  found  from  a  z- velocity  input.  We  relate  the  acoustic  velocity  and  acoustic  potential 
with  the  following  relationship; 

v(x,y,z,t)     =     -V<t>(x,y,z,t).  (2.7) 

Hence,  the  z- velocity  component  is  given  by 

d<p{x,y,z,t) 

v2(x,y,z,t)     = .  (2.8) 

oz 

Since  the  wave  equation  is  in  terms  of  c  (the  acoustic  velocity  in  the  media)  and  time  t,  the 
partial  derivative  with  respect  to  z  must  be  related  to  these  two  parameters.  This  is  done 
by  using  the  fact  that  a  wave  traveling  in  the  positive- z  direction  has  an  argument  of  the 
form  ct  -  z,  resulting  in  the  relationship 

d  d 

Tz  =  ~cTf  (2'9) 

Applying  Eq.  2.9  to  Eq.  2.8  gives  the  z-velocity  at  the  input  plane  (z  =  0)  as 

d<t>(x,y,0,t) 
v2(x,y,z,t)    =     c — .  (2.10) 


Equation  2.10  requires  the  acoustic  potential  at  the  input  plane.  It  was  assumed  in 
Eq.  2.1  that  the  2- velocity  is  separable  which  also  implies  a  separable  acoustic  potential. 
Such  an  acoustic  potential  has  the  general  form 

<f>{x,y,Q,t)  =  s{x,y)r{t)  (2.11) 

at  the  z  —  0  plane.  If  Eq.  2.11  is  substituted  into  Eq.  2.10  and  the  partial  derivative  is 
carried  out,  then  Eq.  2.10  becomes 

v2(x,y,z,t)  =  cs(x,y)r'(t),  (2.12) 

where  the  prime  indicates  a  derivative  with  respect  to  the  time  variable  t.  A  comparison 
of  the  z- velocities  given  in  Eqs.  2.1  and  2.12  indicates  that  T(t)  is  equivalent  to  CT'(t). 
Equation  2.12  is  now  the  input  in  Fig.  2.2c. 

As  stated  earlier,  the  general  solution  of  Fig.  2.2c  is  the  Green's  function.  For  the 
standard  wave  equation  (for  lossless  media)  and  rigidly  baffled  boundary  conditions,  the 
applicable  Green's  function  is  [4] 

h(x,y,z,t)    =     6{C*~RZ)  (2-13) 

where  R  =  \Jx2  +  y2  +  z2.  Substituting  this  into  Eq.  2.2  provides  the  spatial  impulse 
response 


p{x,y,z,t)     =     s{x,y)mx'yh(x,y,z,t) 

6(ct-  z) 
2itR 


=     s(x,y)-x'y    \_„   '  .  (2.14; 


Thus,  the  spatial  impulse  response  in  the  spatial  transfer  domain  is  the  product  of 
the  angular  spectrum  of  the  source  s  and  the  propagation  transfer  function  p,  written 
symbobcally  as 

P{fxJy,Z,t)     =      s{fxJy,t)h{fxJyiZ,t).  (2.15) 

Taking  the  two-dimensional  spatial  Fourier  transform  of  the  Green's  function  in  Eq.  2.13 
gives  the  propagation  transfer  function 

=    2J0  (p\/c2t2  -  z2)  H(ct  -  z) ,  (2.16) 

where  the  relationship 


6{t-\)\  Jo  {p^/cH2  -  z2)  H(ct  -  z) 

Rn      J     "  (ct)n~l 


(2.17) 


was  applied  (with  r  =  .//J  +  /*).  The  term  H (ct  —  z)  is  the  Heaviside  step  function  and 
makes  the  wave  causal. 

Equation  2.16  exhibits  two  important  points.  The  first  is  that  the  propagation  transfer 
function  is  a  Bessel  of  the  first  kind  of  zero  order.  Secondly,  the  propagation  transfer 
function  can  be  identified  as  a  time-varying  spatial  filter  having  a  Bessel  shape.  The  filter 
begins  as  an  all-pass  filter  early  in  its  time  evolution  and  then  becomes  more  of  a  low-pass 
filter  as  time  progresses.  This  increasingly  narrow  low-pass  filter  reduces  the  resolution 
in  the  receiving  plane  of  the  system  as  time  progresses.  In  our  computer  programs,  the 
propagation  spatial  filter  at  various  instants  of  time  are  produced  by  the  MATLAB  file, 
ACJTL.M. 

Equation  2.16  is  valid  for  an  input  that  is  temporally  impulsive  and  spatially  arbitrary. 
Taking  the  inverse  two-dimensional  spatial  transform  of  Eq.  2.16  would,  in  this  case,  give 
a  final  result  since  convolution  with  an  impulse  is  the  same  function  at  the  time  that  the 
impulse  occurred.  To  account  for  a  non-impulsive  time  input  component,  the  convolution 
of  Eq.  2.5  must  be  carried  out,  resulting  in 


<f>(x,y,z,t)     =    T{tytT-'{f>UxJy,z,t))  ,  (2-18) 

where  p  is  the  product  of  the  angular  spectrum  of  the  source  5  and  the  propagation  trans- 
fer function  h.  Since  only  temporally  impulsive  inputs  are  simulated  in  the  examples  that 
follow,  the  convolution  in  Eq.  2.18  was  not  computed  for  our  cases.  Our  final  solution, 
therefore,  reduces  to  taking  the  inverse  2-D  spatial  Fourier  transform  of  the  spatial  im- 
pulse response.  The  program  module,  AC-PROP. M,  simulates  the  spatial  impulse  response 
solution  for  the  input  excitation  function  distribution  chosen  by  the  user. 

The  program  that  implements  these  equations  is  discussed  in  Chapter  III  and  detailed 
in  Appendix  C.  Illustrative  examples  of  the  time-varying  filters  and  outputs  follow  in  the 
Chapter  IV  discussion  with  more  examples  supplied  in  Appendices  B  and  E.  The  tool 
of  implementation  for  the  program  was  MATLAB  with  graphical  assistance  provided  by 
AXUM.  Output  data  was  also  represented  using  a  scientific  visualization  program  called 
Spyglass  Dicer.  An  overview  of  these  programs  follows  in  the  next  section. 

2.4      Software  Tools 

2.4.1      MATLAB  Overview 

MATLAB  is  a  high-performance,  interactive  numeric  computation  software  package  for 
science  and  engineering  applications  produced  by  The  MATH  WORKS,  Inc.  The  name 
comes  from  MATrix  LABoratory;  hence,  the  basic  data  element  is  a  matrix  (or  array), 
which  does  not  require  dimensioning.  MATLAB.  A  major  advantage  of  MATLAB  is  that 
it  uses  a  "vectorized"  approach  to  computations,  simplifying  the  programming.  Another 
distinct  advantage  is  the  availability  of  preprogrammed  functions,  such  as  calculation  of  the 
two-dimensional  FFTs  and  the  Bessel  function  [16]. 


In  MATLAB  there  are  two  type  of  macro-like  files.  A  script  file  is  used  to  automate  long 
sequences  of  commands  including  functions.  Arguments  are  not  passed  into  script  files.  A 
function  file,  however,  may  have  arguments  passed  into  it  and  out  of  it.  Examples  of  script 
files  in  this  thesis  include  AC-FIL.M  and  AC-PROP.M.  Examples  of  functions  include  the 
input  functions,  the  three  dimensional  graphing  function  mesh,  and  the  two  "fit"  functions 
that  realize  the  Fourier  transform  [16]. 

The  two  "fft"  functions  employed  for  the  Fourier  transform  are  fft2  and  fftshift.  The 
Aft 2  function  is  a  two-dimensional  fast  Fourier  transform  and  is  repeatedly  used  throughout 
our  program.  MATLAB  allows  the  use  of  only  positive  indices  in  an  array.  Geometrically, 

this  is  equivalent  to  working  only  in  the  first  quadrant  of  an  x y  plane.  Our  sources  are 

symmetric  around  the  origin  and  we  want  to  display  them  using  all  four  quadrants.  The 
key  to  the  transformation  is  the  fftshift  function.  Figure  2.3a.  shows  the  first-quadrant 
representation  of  a  rectangular  source  centered  on  the  origin.  (Note  the  MATLAB  assumes 
implicitly  that  the  function  is  periodic  in  the  x  and  y  directions.  Periodic  translation  of  the 
first-quadrant  in  the  x  and  y  directions  will  fill  in  the  rest  of  the  source  at  the  origin.)  The 
fftshift  function  rearranges  the  data  so  that  it  is  centered  in  the  first  quadrant.  (Actually 
the  data  is  not  quite  centered  when  one  uses  an  even  number  of  sample  points,  since  there  are 
not  any  sample  points  aligned  with  the  center.  The  shifted  array  is  offset  from  the  quadrant 
center  by  one-half  of  a  sample  interval  in  both  the  x  and  y  directions.)  Figure  2.3b  represents 
the  shifted  source.  For  viewing  purposes,  the  shifted  version  is  easier  to  understand;  for 
computational  purposes  the  unshifted  function  must  be  used  to  avoid  phase  errors  when 
calculating  the  transform  or  its  inverse.  We  will  use  the  terms  "shifted"  for  functions 
represented  in  the  centered  geometry  and  "unshifted"  for  the  functions  represented  in  the 
corner  geometry.  The  interested  reader  is  referred  the  MATLAB  User's  Guide  [16]  for  more 
details. 

The  MATLAB  mesh  command  can  be  used  to  provide  three-dimensional  plots  of  the 
computed  fields.  Alternatively,  the  input  and  output  data  can  be  stored  in  ASCII  format, 
imported  into  a  more  powerful  graphics  package,  called  AXUM,  and  plotted. 

2.4.2     AXUM  Overview 

AXUM  is  an  interactive  software  package  for  technical  graphics  and  data  analysis,  produced 
by  TriMetrix  Inc.  AXUM  allows  easy  manipulation  of  raw  data  imported  from  MATLAB 
in  ASCII  format.  (Data  may  also  be  imported  from  several  other  formats.)  Once  imported, 
data  manipulation  can  be  carried  out  by  AXUM  through  use  of  the  Transform  and  Convert 
menu  options.  Then  the  data  is  arranged  using  the  MAT2GR1D  routine  found  in  the  History 
Editor  menu,  so  that  the  desired  surface  plot  can  be  generated.  The  MAT2GR1D  routine 
produces  three  columns  of  data  from  the  original  data  array  by  placing  each  element  of 
acoustic  potential  data  in  the  z  data  column  with  the  indices  from  the  array  in  the  x  and 
y  columns.  Once  processed,  the  data  is  ready  for  graphing  [17]. 

The  Graph  menu  gives  various  options  for  controlling  the  graph's  attributes  and  general 
aesthetics.  Once  the  desired  axes  intervals,  labels,  and  titles  have  been  set,  the  graph  can 
be  displayed  in  the  Edit  Screen  window.    The  Edit  Screen  window  allows  changes  to  be 
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made  to  the  graph  while  it  is  displayed.  (This  can  be  is  time-consuming  because  the  graph 
is  redrawn  after  each  change.)  A  useful  feature  of  the  Edit  Screen  window,  however,  is  the 
capability  to  rotate  and  tilt  the  graph  interactively  so  that  the  preferred  aspect  is  achieved. 
After  establishing  the  desired  attributes,  the  graph  can  be  saved  as  a  graph  or  as  an  image. 
The  concept  of  saving  the  graphs  as  images  is  another  very  useful  feature  because  it  allows 
a  single  graph  to  be  assembled  and  saved  as  a  graph  template  on  which  the  other  images 
may  be  overlaid. 

2.4.3     Spyglass  Software  Overview 

The  Spyglass  software  consists  of  a  set  of  commercial  programs  produced  by  Spyglass,  Inc. 
for  the  purpose  of  aiding  in  visualization  of  scientific  data  on  Macintosh  computers.  (One 
program,  the  Transform  program,  is  also  commercially  available  for  Unix  workstations.) 

The  Spyglass  Data  Utility  has  the  capability  to  read  data  from  a  wide  variety  of  formats 
into  the  HDF  (Hierarchical  Data  Format)  format  1  used  by  the  Spyglass  software.  In 
particular,  for  our  work,  the  utility  is  able  to  read  ASCII  data.  In  our  case,  we  took  the  64 
ASCII  data  files  produced  by  propagation  model  and  transferred  them  to  a  Macintosh  II 
computer.  The  Spyglass  Data  Utility  was  used  to  create  an  HDF  file,  representing  a  data 
volume  of  128x128x64  samples.  (The  typical  data  compression  was  from  16  MB  of  ASCII 
data  to  a  6-MB  HDF  file.) 

Spyglass  Dicer  is  a  program  that  allows  the  user  to  interactively  look  at  the  3-D  data 
volume  in  various  ways.  Planes  of  any  arbitrary  orientation  can  be  inserted  or  rectan- 
gular regions  of  the  data  can  be  observed.  (Figures  4.3  and  4.4  on  pages  24  and  25  are 
representative  of  the  data  views  that  can  be  obtained  with  the  Dicer  program.) 

The  Spyglass  Transform  program  allows  the  user  to  pick  a  plane  of  the  HDF  data  volume 
(currently,  the  planes  must  be  oriented  perpendicular  to  one  of  the  axes)  and  to  view  it  in 
a  variety  of  formats,  such  as  a  color  intensity  display,  a  gray-scale  intensity  display,  or  a 
surface  plot. 


JThe  HDF  format  is  becoming  a  popular  means  of  transferring  large  data  files  between  computinj 
platforms. 
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Chapter  3 

MATLAB  Modeling  of  Equations 


This  chapter  discusses  the  two  MATLAB  files  that  implement  the  model,  AC-F1L.M  and 
AC-PROP.M.  Together,  these  script  files  form  the  two  modules  of  the  program  that  simulate 
acoustic  wave  diffraction  via  implementation  of  the  concepts  and  equations  presented  in  the 
previous  chapter.  The  program  was  written  in  two  modules  since  the  propagation  filter 
characteristics  of  the  medium  depend  only  on  the  locations  of  the  source  and  receiving 
plane  locations  and  are  independent  of  the  source  geometry.  The  propagation  filters  could 
be  computed  once  and  then  be  reused  for  a  variety  of  excitation  conditions.  This  modularity 
of  approach  was  fortuitous,  since  the  calculation  of  the  Bessel  functions  in  the  filters  was 
computationally  intensive.  (This  modular  approach  to  the  calculation  of  the  fields  is  one 
of  the  major  advantages  of  our  method.)  A  working  narrative  of  AC-FIL.M  opens  the 
discussion,  followed  by  a  working  narrative  of  AC-PROP.M.  The  excitation  functions,  or 
input  functions,  are  included  under  the  AC-PROP.M  heading.  A  brief  summary  of  the 
program  steps  follows  in  the  final  section. 

3.1     Acoustic  Filter  Module 

The  script  file,  AC_FIL.M  (ACoustic  FILter)  computes  the  time-varying  Bessel  filters  of 
Eq.  2.16,  repeated  here  for  convenience  as 

=    2 Jo  (pVc2t2  -  z2\  H(ct -  z) .  (3.1) 

Before  discussing  the  coding  of  Eq.  3.1,  the  array  geometry  must  be  explained. 

The  array  geometry  is  set  by  the  size  of  the  array  and  the  sampling  frequency.  The 
variable  base  denotes  the  number  of  points  on  a  side  in  the  base  array  creating  an  base  x  base 
array.  (It  is  advantageous  in  calculating  an  FFT  to  make  base  a  power  of  2;  this  is  not 
necessary,  however.)  Initially  base  was  given  a  value  of  64  as  in  previous  work  [3] —  [7]. 
Once  the  program  result  was  verified,  base  was  increased  by  a  factor  of  two  to  128.  Making 


12 


base 


bose 


+  1 


Suborroy  t 

Suborroy  1 

Suborray  HI 

^  (N0.N0) 
Suborroy  IV 

base 


base 


Array  index 
Figure  3.1:  Offset  geometry  of  base  array  matrix. 


base  an  even  number,  however,  causes  the  center  of  symmetry  of  the  array  to  fall  between 
array  elements.  The  center  of  symmetry  of  the  source  A70  coincided  with  the  array  element 
at 


base 

NO    =     — —  +  1. 


(3.2) 


This  results  in  an  offset  geometry  as  shown  in  Figure  3.1. 

The  offset  center  at  (NO,  NO)  divides  the  base  array  into  four  subarrays,  each  of  different 
sizes.  The  fact  that  the  subarrays  have  different  sizes  is  important  in  the  use  of  symmetry 
because  only  one  subarray  is  actually  computed;  the  other  subarrays  are  determined  by 
symmetry.  We  calculate  the  data  for  the  NO  x  (NO  -  1)  points  in  Subarray  I,  shown  in 
Fig.  3.1.  In  addition,  data  is  calculated  for  an  additional  column  (WO  x  1)  located  at  the 
right  side  of  Subarray  I.  (The  extra  column  was  required  to  compute  all  of  the  values  in  the 
other  three  quadrants  using  symmetry.)  Subarrays  II,  III,  and  IV  of  Fig.  3.1  are  determined 
from  Quadrant  I  using  symmetry  with  MATLAB's  flip  commands.  The  use  of  symmetry  in 
this  manner  was  employed  in  both  program  modules.  It  should  be  noted  that  this  assumed 
x-  and  j/-axis  symmetry  imposes  an  other  restriction  on  the  assumed  source  geometry;  this 
restriction  can  be  removed  by  performing  the  computations  over  the  full  dimensions  of  the 
array  rather  than  relying  on  symmetry  to  simplify  the  calculations. 

The  number  of  time  samples  (or  number  of  time  slices),  slices,  was  set  at  64.  Initially, 
this  was  to  emulate  previous  work  [3] —  [7];  however,  the  time  duration  and  resolution 
proved  to  be  adequate  for  the  follow-on  simulations.    Although  there  are  64  time  slices, 
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only  61  filters  are  generated  by  the  MATLAB  implementation  of  Eq.  3.1.  The  Heaviside 
step  function  in  Eq.  3.1  was  simulated  by  arbitrarily  setting  a  variable  Step  to  three  which 
produces  three  all-zero  rows  for  the  first  three  time  slices.  The  result  is  that  slices  —  Step, 
or  61,  time  slices  actually  require  computing.  Since  the  first  three  time  slices  are  zero,  the 
computations  start  with  the  fourth  time  slice  at  time  z/c  and  proceed  to  the  time  value, 
time.max,  provided  by  the  user. 

Referring  back  to  Eq.  3.1,  it  is  seen  that  the  argument  to  the  Bessel  function,  p\Zc7t2  —  z2, 
is  composed  of  four  variables.  Of  the  four,  only  the  time  t  varies  within  the  program.  In 
the  MATLAB  code,  time  t  is  represented  by  the  variable  time.  As  previously  stated,  filter 
generation  does  not  start  until  time  z/c  and  is  linearly  incremented  to  the  maximum  time 
of  propagation  time.max.  The  source-to-receiver  distance  z  was  assigned  the  value  2  =  10 
cm  in  the  MATLAB  code  and  c,  the  acoustic  velocity  in  the  medium,  was  assigned  the  value 
c  =  1500  m/s  (velocity  in  water).  The  value  of  z  was  originally  chosen  to  parallel  the  pre- 
vious work  [3] —  [7]  and  was  found  to  be  convenient  for  Subsequent  simulations.  The. value 
of  timejmax  was  set  to  150  ps  for  the  verification  phase  and  to  375  ps  for  the  remainder 
of  the  simulations.  Consequently,  time  ranged  from  66.667  ps  to  150  ps  (or  375  ps)  in  61 
increments. 

The  final  variable  in  the  Bessel  argument  to  be  examined  is  p.  In  the  MATLAB  code,  p 
was  given  the  name  rho  and  has  a  maximum  value,  rho.max,  of  200.  The  value  of  rho.max 
=  200  was  arbitrarily  chosen  subject  to  the  constraint  that  it  needs  to  be  large  enough  that 
the  field  does  not  extend  beyond  the  edges  of  the  array  at  its  largest  width  (or  else  the  field 
will  be  aliased  back  into  the  array  from  the  edges).  The  value  of  200  was  chosen  following 
Merrill  [7].  Although  rho  is  not  time-varying,  it  does  vary  with  spatial  frequency 

P    =     v//x2  +  /y2-  (3-3) 

A  vector  having  NO  —  1  points  extending  from  0  to  rho.max  was  then  formed.  Then, 
the  vector  was  used  in  the  MATLAB  routine  meshdom  [16]  to  form  two  identical  matrices 
rhox  and  rhoy.  The  meshdom  routine  takes  a  given  vector  and  forms  two  matrices  with  the 
same  spacing  in  the  x  direction  and  y  direction.  The  two  matrices,  rhox  and  rhoy  represent 
fx  and  fy  in  Eq.  3.3  and  Fig.  3.2,  which  shows  graphically  the  construction  of  rho.  In 
Fig.  3.2,  the  inner  scales  are  in  terms  of  the  column  and  row  numbers,  respectively.  The 
row  index  runs  from  1  to  NO;  the  column  index  runs  from  NO  to  base.  The  x  and  y  axes 
can  be  rescaled  to  units  of  spatial  frequency  (fx  and  fy)  with  units  of  cycles/m.  This  is 
represented  in  Fig.  3.2  by  the  outer  labeling. 

The  combination  of  rho  and  the  other  variables  forms  the  argument  arg  to  the  MATLAB 
Bessel  function.  Since  time  varies  for  each  time  slice,  arg  varies  for  each  time  slice  generating 
a  filter  per  time  slice.  After  generation  each  filter  is  stored  to  disk  for  use  by  ACPROP.M. 
The  variables  base,  NO,  slices,  and  Step  are  also  stored  to  disk  for  use  by  ACPROP.M. 
The  interested  reader  can  find  a  detailed  explanation  and  the  source  code  in  Appendix  A. 
Graphical  examples  of  the  time- varying  filters  are  include  as  Appendix  B. 
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Figure  3.2:  Construction  of  rho  shown  graphically. 


3.2     Acoustic  Propagation  Module 

The  script  file,  AC_PR0P.M  (ACoustic  PROPagation),  allows  the  user  to  specify  the  type 
and  dimensions  of  the  input  excitation.  It  then  takes  the  user's  chosen  input  function 
and  simulates  propagation  as  a  function  of  time.  AC-PROP  computes  the  spatial  impulse 
response  p(x,  y,  z,  t), 


p{x,y,z,t)     =    ^/v{s(/x,/v)/i(/x,/j,,2,*)}  , 


(3.4; 


using  the  time-varying  Bessel  filters,  produced  by  ACJFIL  for  h(fx,fy,z,t).  In  Eq.  3.4. 
Tl  f  {•}  is  the  inverse  spatial  Fourier  transform  operator  and  the  denotes  a  transformed 
function. 

To  compute  p(x,y,z,t),  AC_PR0P  first  loads  the  variables  passed  from  AC.FIL  and 
then  queries  the  user  to  make  an  input  function  choice.  The  user  is  given  four  input 
functions  from  which  to  choose:  Circle,  Table,  Gaussian,  and  Bessel.  The  Table  and  Circle 
are  equal-amplitude  sources  having  the  shape  of  a  square  and  a  circle,  respectively  (known 
in  acoustics  as  the  square  and  circular  piston).  The  Gaussian  and  Bessel  inputs  are  circularly 
truncated  functions  that  have  the  indicated  amplitude  distribution  within  the  circle  (i.e., 
the  amplitudes  are  spatially  varying  across  the  face  of  the  source).  Pursuant  to  the  input 
function  choice,  the  user  is  asked  to  input  the  diameter  d  of  the  truncating  circle  (or  width  of 
the  truncating  square  w  in  the  case  of  the  Table  function).  Once  the  program  was  verified, 
a  diameter  d  of  51  was  used;  this  value  is  available  as  the  default  in  the  menus.    In  the 
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case  of  the  Gaussian  or  Bessel  input,  the  user  is  further  requested  to  input  the  standard 
deviation  sigma  or  a  scaling  factor  a,  respectively.  These  two  function  defining  parameters 
were  varied  on  a  case-by-case  basis. 

The  variable  shft.input  holds  the  computed  input  array  that  was  generated  by  the  func- 
tion written  to  model  the  chosen  excitation  function  (Appendix  D).  (The  reader  is  reminded 
that  shft  prefix  in  a  variable  name  indicates  the  array  has  the  centered  geometry  as  dis- 
cussed in  Chapter  II.)  From  shft.input,  F.input  is  created  by  shifting  (fftshift)  shft.input 
to  a  corner  geometry  and  then  taking  the  two-dimensional  spatial  Fourier  transform  (fft2). 
As  explained  in  Chapter  II,  the  fftshift  operation  is  necessary  before  the  Fourier  trans- 
form operation  to  obtain  the  correct  phase  relationship  in  the  transform  operation.  With 
a  shift  back  to  the  center  geometry,  the  angular  spectrum  of  the  source  s(/x,/y),  called 
shft.F.input  in  the  program,  is  created  for  user  viewing,  if  desired.  The  Bessel-function 
propagation  transfer  function  h(fx,  fy,z,t)  from  Eq.  3.1  must  now  be  loaded  so  that  the 
angular  spectrum-propagation  transfer  function  product  sh  on  the  right  side  of  Eq.  3.4 
can  be  formed.  The  loading  and  multiplication  process  is  repetitive  since  shft.F.input  must 
form  a  product  with  the  filter  (or  appropriate  propagation  transfer  function)  from  each  time 
slice.  This  repetitive  multiplication  is  accomplished  with  a  loop. 

The  product  of  s  and  h  (called  shft. F. output).  To  find  the  desired  result,  the  two- 
dimensional  inverse  spatial  Fourier  transform  ( ifft 2 )  must  be  taken.  Before  this  can  be 
done,  F.output  is  formed  by  shifting  shft.F.output  from  the  centered  geometry  to  the  corner 
geometry.  Executing  the  inverse  transform  of  the  product  yields  the  output  (output)  which 
is  then  shifted  to  give  shft. output  for  viewing.  The  array  shft.output  represents  the  output  at 
the  time  slice  that  the  loop  is  currently  computing;  shft.output  does  not  depict  the  acoustic 
potential  (or  propagation  pattern)  through  time;  it  only  depicts  the  acoustic  potential  at  a 
specific  time. 

To  produce  a  time  history  of  the  desired  output,  the  center  row  (row  NO)  of  shft.output 
is  taken  and  placed  in  the  array  output.plot  as  the  m-th  column  (m  is  the  loop  counter 
which  relates  directly  to  the  time  slice  number,  i.e.,  when  m  =  4,  the  computations  for  the 
fourth  time  slice  are  performed).  This  results  in  an  array  whose  size  is  base  x  slices  when 
Step  zero-valued  rows  are  added  (that  precede  time  slice  NO  row.  The  output  examples 
are  graphical  interpretations  of  output.plot.  Results  generated  in  this  manner  are  given 
in  Chapter  rV  for  all  of  the  excitation  functions.  A  detailed  explanation  of  AC-PROP  is 
provided  in  Appendix  C. 

3.3      Program  Summary 

The  previous  two  sections  gave  an  overview  of  the  two  program  modules,  AC-FIL  and 
AC_PROP,  including  the  code  variable  names  and  values  assigned.  What  follows  here  is  a 
summary  of  the  steps  that  the  program  accomplishes.  Step  one  is  accomplished  by  AC-FIL. 
In  this  step  the  slices-Step  filters  to  be  used  by  AC-PROP  are  generated  and  saved. 

AC-PROP  generates  the  user-specified  excitation  function  s(x,y).  Then  the  angular 
spectrum  of  the  source  s(fx,fy)  is  computed  by  taking  the  2-D  spatial  Fourier  transform 
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of  s(x,y).  The  product  sh  is  computed  for  each  time  slice  via  a  loop.  In  the  loop,  the 
inverse  2-D  transform  of  the  product  forms  an  output  for  the  specific  time  slice.  The  JVO-th 
(center)  row  of  that  output  is  then  placed  in  successive  columns  of  a  new  array  to  form  the 
final  output,  p(x, 0, 10,  t). 

In  the  following  chapter,  examples  of  the  final  outputs  are  given.  The  excitations  used 
for  verification,  as  previously  related,  are  the  Table  and  Circle  excitation  functions.  New 
values,  as  explain  in  this  chapter,  were  then  used  for  the  simulation  of  the  Gaussian  and 
Bessel  excitations. 
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Chapter  4 

Numerical  Simulations 


In  the  previous  chapter,  a  functional  explanation  of  the  two  program  modules  was  given 
including  values  assigned.  The  first  section  of  this  chapter  reiterates  the  defining  parameters, 
gives  a  brief  explanation  of  each,  and  gives  the  parameter's  units.  In  the  following  section, 
the  defining  parameters  are  given  the  values  used  to  verify  correct  operation  of  the  program 
for  the  circular  and  square  piston  sources.  The  last  section  presents  results  for  the  non- 
piston  Gaussian  and  the  Bessel  excitation  functions. 

4.1      Defining  Parameters 

A  defining  parameter  is  a  parameter  that  delineates  an  aspect  of  the  basic  setup  upon 
which  all  the  remaining  parameters  or  variables  depend.  In  the  work  of  this  thesis  there 
are  two  sets  of  defining  parameters — those  for  the  filter  generation  and  those  for  generation 
of  the  excitation  functions.  The  filter  parameters  are  found  at  the  beginning  of  AC-FIL. 
AC-PROP  reads  these  parameters  from  a  data  file  that  is  stored  with  the  results  of  the 
filter  calculations. 

The  defining  parameters  found  in  ACJPIL  include  base,  slices,  Step,  c,  z,  time.max,  and 
rho.max.  The  first  parameter  base  sets  the  dimensions  of  the  base  array,  giving  the  number 
of  sample  points.  The  dimensions  of  the  base  array  are,  therefore,  base  x  base  where  base  is 
required  to  be  a  power  of  two  (typically,  128).  Making  base  a  power  of  two  allows  MATLAB 
to  use  a  high-speed  radix-2  fast  Fourier  transform  algorithm  [16]  to  compute  the  spatial 
transforms.  The  next  parameter  slices  is  the  number  of  time  samples.  Of  these  slices  time 
slices,  slices-Step  slices  require  filters  to  be  computed;  the  parameter  Step  is  the  number  of 
leading-zero  rows  in  the  base  x  slices  output  array;  as  explained  in  the  preceding  chapter, 
this  simulates  the  Heaviside  step  function.  The  parameters  base,  slices,  and  Step  are  unitless 
and  are  stored  by  AC-FIL  in  a  file  for  use  by  AC-PROP. 

The  remaining  defining  parameters  of  AC-FIL  have  units  and  are  used  only  in  the 
computations  of  the  filters.  The  acoustic  velocity  in  the  medium,  free-space  in  this  case, 
is  denoted  by  the  parameter  c  having  the  units  of  m-s-1.  The  source-to-reception  point 
distance  has  the  designation  z  with  units  of  meters.    The  maximum  time  of  propagation 
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Parameter 

Definition  (units) 

base' 

Size  of  square  base  array 

slices' 

Number  of  time  slices 

Step' 

Number  of  leading  zero-rows 

c' 

Acoustic  velocity  in  media  (m/s) 

f 

Distance,  source-to-receiver  (m) 

time.max" 

Maximum  time  of  propagation  (s) 

rho.max' 

Spatial  radius  of  the  filters  (length-1) 

'  passed  to  AC_PROP.M 

Table  4.1:  Summary  of  the  defining  parameters  in  AC-FIL. 


time.max  has  units  of  seconds.  The  spatial-frequency  radius  of  the  filters  rho.max  (or  rho) 
has  units  of  inverse  length  (i.e.,  m-1,  cm-1,  etc.).  The  unit  of  length  depends  on  the  area  to 
be  represented  by  the  base  array.  These  four  parameters  relate  directly  to  Eq.  3.4  and  are 
the  parameters  that  dictate  the  diffraction  properties  of  the  filters.  Table  4.1  summarizes 
the  defining  parameters  found  in  ACJPIL. 

Another  important  set  of  defining  parameters  is  the  set  that  defines  the  user  chosen  in- 
put. These  input  defining  parameters  are  entered  by  the  user  when  requested  by  AC_PROP. 
Once  the  input  function  is  chosen,  the  diameter  of  the  truncation  circle  d  is  input.  (The 
width  of  the  table  w  (vice  d)  is  input  in  the  case  of  the  Table  excitation.)  The  parameters 
d  (or  w)  are  expressed  as  the  number  of  points,  out  of  base  total  points,  that  define  the 
diameter  (or  width)  of  the  function 

To  transition  from  a  diameter  in  terms  of  a  number  of  points  to  an  actual  metric  value, 
two  equations  were  needed.  The  equations  are 


Ax    = 


1 


2Pr 


and 


d    = 


k&x 
k 

•^Pmaj 


(4.1) 

(4.2) 
(4.3) 


where  Ax  is  the  length  of  a  segment,  pm^x  is  the  maximum  spatial  radius,  k  is  the  number 
of  segments,  and  d  is  the  diameter  (or  width  w  for  the  Table  function).  To  determine  the 
metric  diameter,  Eq.  4.1  was  used  to  solve  for  Ax  by  setting  pmAX  to  200  m-1.  This  value 
resulted  in  a  Ax  of  2.5  x  10-3  m  or  2.5  mm. 

If  the  Gaussian  excitation  is  selected,  the  user  enters  the  value  of  the  standard  deviation 
a,  upon  request.  In  a  Bessel  excitation  selection  the  user  enters  the  scaling  factor  a  when 
requested.  Table  4.2  gives  a  summary  of  the  defining  parameters  used  in  AC_PROP. 
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Parameter 

Definition  (units) 

base" 

Size  of  square  base  array 

slices" 

Number  of  time  slices 

Step' 

Number  of  leading  zero-rows 

c' 

Acoustic  velocity  in  media  (m/s) 

? 

Distance,  source-to-receiver  (m) 

time.max' 

Maximum  time  of  propagation  (s) 

rho.max" 

Maximum  spatial  frequency  (m-1) 

d 

Diameter  of  excitation  function  (sample  points) 

w 

Width  of  Table  excitation  function  (sample  points) 

a 

Gaussian  standard  deviation 

a 

Bessel  scaling  factor 

*  passed  from  AC-FIL.M 

Table  4.2:  Summary  of  the  defining  parameters  used  in  AC-PROP. 


4.2      Program  Verification 

In  verifying  the  program  output,  two  excitation  functions  were  used,  the  Table  and  the  Circle 
functions.  The  outputs  generated  by  the  program  from  these  excitations  were  compared 
to  the  results  found  with  the  previous  FORTRAN  implementations  [4,  3,  7]  for  validation. 
After  a  general  explanation  about  the  generation,  formatting,  and  titling  of  the  outputs, 
the  two  excitation  functions  are  presented.  The  Table,  the  first  verification  function,  is  then 
discussed  and  a  table  of  defining  parameters  is  given.  The  second  verification  function,  the 
Circle,  follows  with  a  similar  discussion  and  table  of  defining  parameters. 

4.2.1     Format  of  Results 

The  graphical  outputs  for  the  two  excitation  functions  used  for  verification,  the  Table  and 
the  Circle,  were  generated,  formatted,  and  titled  in  the  same  manner.  The  outputs  are  for 
a  source-to-receiver  distance  of  z  =  10  cm.  Each  output  was  for  a  given  time  slice  and 
consisted  of  a  128x128  array  of  values.  There  were  64  time  slices  including  3  leading  all-zero 
arrays  which  simulate  the  step  function  at  the  arrival  of  the  wave  (t  =  z/c).  The  MATLAB 
programs  contain  optional  commands  to  plot  the  output  data  using  MATLAB 's  plotting 
routines  or  to  store  the  data  in  ASCII  format  for  importation  into  the  AXUM  plotting 
program.  (The  AXUM  program  is  more  flexible  than  the  MATLAB  plotting  routines.) 

Additionally,  the  64  data  arrays  were  combined  with  the  Spyglass  Data  Utility  into  a 
128x128x64  array.  Spyglass  Transform  and  Dicer  could  be  used  to  interactively  pick  the 
data  slice  of  interest  and  to  plot  the  results. 
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NAME 

VALUE 

DEFINITION 

base 

128 

Size  of  square  base  array 

slices 

64 

Number  of  time  slices 

Step 

3 

Number  of  leading  zero-rows 

c 

1500 

Acoustic  velocity  in  media  (m/s) 

z 

0.1 

Distance,  source-to-receiver  (m) 

time.max 

3.75  x  10~4 

Maximum  time  of  propagation  (s) 

rho.max 

200 

Spatial  radius  of  the  filters  (length-1) 

w 

23,31 

Width  of  Table  (samples) 

Table  4.3:    Values  of  defining  parameters  for  the  Table  input  function  used  for  program 
verification. 


4.2.2      Table  Impulse  Excitation 

The  first  excitation  function  to  be  run  by  the  program  was  the  Table  function  (for  a  square 
piston  source).  There  were  several  reasons  to  use  the  Table  as  the  first  input;  the  Table 
function  is  an  easy  function  to  implement  and  the  results  could  be  readily  compared  to 
results  found  in  the  literature  [3,  5,  6,  7].  Table  4.3  provides  list  of  the  defining  parameters 
used. 

The  values  of  base,  slices,  z,  time.max,  rho.max,  and  w  chosen  in  Table  4.3  parallel 
those  found  in  the  literature  used  for  validation.  The  acoustic  velocity  c  of  1500  m/s  is 
the  velocity  in  water.  To  simulate  the  step  function,  a  number  Step  of  leading  zero-rows 
was  incorporated  and  arbitrarily  set  to  three.  The  results  were  comparable  to  those  in  the 
literature  and  are  presented  (with  the  input  functions)  in  Figs.  4.1  and  4.2. 

Here  the  widths  w  —  23  and  w  =  31  samples  translate  into  metric  widths  of  w  —  5.5  cm 
and  w  —  7.5  cm,  respectively.  In  the  w  =  23  case,  k  =  22  was  the  input  in  Eq.  4.3  and,  in 
the  w  =  31  case,  the  input  was  k  =  30  samples.  Note  that  the  x  and  y  axes  of  the  input 
functions  range  from  1  to  128,  delineating  a  128x128  base  array. 

The  outputs  present  the  magnitude  of  the  acoustic  potential  at  an  observation  point  10 
cm  from  the  source,  p(x,0, 10,  t).  A  diffraction  duration  from  the  initial  impulse  (t  =  0)  to 
t  —time.max  (150  /xs)  is  represented  as  a  function  of  radial  distance;  i.e.,  p(x, 0,10,0)  to 
p(x, 0,10,150)  is  represented.  This  gives  the  3-D  view  of  the  general  diffraction  through 
time.  The  output  images  show  several  interesting  features. 

The  development  of  "tails,"  explained  in  terms  of  edge  waves  [3,  18],  can  be  seen  in 
the  3-D  images.  Also  of  note  are  the  overshoots,  having  a  maximum  magnitude  of  2.11 
(both  cases),  and  the  undershoots  (difficult  to  see  in  these  two  cases);  these  are  due  to  the 
additive  nature  of  the  interference  patterns  of  the  waves  originating  on  the  edges  of  the 
discontinuous  source. 
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(o) 


(b) 


Figure  4.1:  Table  spatial  input  and  time-space  output  for  w  =  23  samples  (w  —  5.5  cm)  at 
z  =  10  cm. 
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(a) 


(b) 

Figure  4.2:   Table  spatial  input  and  time-space  output  for  w  =  31  samples  (w  =  7.5  cm)  at 
z  =  10  cm. 


23 


Propag 


atonof  table(128 


iO 


r     p 

iff 

.: •„  ,„,,.„, .  -^.  «i-  ■  •> ~  s^i*: 

127 


time 


240 


-0.1455 


Amplitude 


2070 


Figure  4.3:  Table  output  for  d  —  31  samples.  Dicer  representation. 


3D  outputs  for  Table  excitation 

The  Spyglass  Dicer  program  allows  selected  three-dimensional  representations  of  the  propa- 
gation output.  Figure  4.3  shows  a  rectangular  cube  set  into  one-fourth  of  the  data  volume. 
The  x  and  y  axes  represent  the  spatial  variables  and  range  from  0  to  127  samples.  The  time 
axis  represents  the  64  time  slices  of  data  that  were  computed.  To  elongate  the  time  axis, 
the  Dicer  program  was  used  to  increase  the  number  of  slices  by  a  factor  of  4x  (to  a  total  of 
256  slices);  linear  interpolation  is  used  to  compute  the  values  of  the  expanded  slices  lying 
between  the  original  slices.  Due  to  the  nature  of  the  interpolation  used,  only  the  expanded 
slices  ranging  from  0  to  252  are  shown  on  the  time  axis.  The  three  visible  surfaces  of  the 
rectangular  cube  show  the  calculated  data  in  three  planes.  The  three  planes  are  located  at 
1  =  64,  y  =  65,  and  time  =  252.  An  additional  plane  is  located  at  time  =  16  to  show  the 
shape  and  size  of  the  source  (since  the  theory  predicts  that  the  field  at  time  =  16  expanded 
samples  [when  r  =  z/c]  duplicates  the  source  excitation). 

An  alternative  Dicer  representation  of  the  field  is  shown  in  Figure  4.4.  Here,  a  horizontal 
slice  of  the  data  is  placed  at  y  =  65  and  five  vertical  slices  of  the  data  are  located  at  x  - 
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Figure  4.4:  Table  output  for  d  =  31  samples.  Alternative  Dicer  representation. 
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16,  64,  132,  200,  and  252.  (The  values  were  arbitrarily  chosen,  but  were  restricted  to  be  a 
multiple  of  four  to  avoid  showing  interpolated  values  of  data.) 

Appendix  E  contains  Dicer  representations  for  the  other  Table  and  Circle  excitations 
that  were  studied. 

4.2.3      Circle  Impulse  Excitation 

The  defining  parameters  for  the  Circle  excitation  are  the  same  as  those  introduced  in  Ta- 
ble 4.3  with  the  exception  that  only  the  d  =  31  case  is  presented.  A  diameter  d  =  31  samples 
translates  into  a  metric  diameter  of  d  =  7.5  cm.  Again  the  results  are  comparable  to  those 
found  in  the  literature  [5,  7].  Figure  4.5  gives  the  input  function  and  the  ensuing  outputs. 
As  in  the  Table  case,  the  base  array  is  a  128x128  array  with  the  propagation  pattern  formed 
by  successive  p(x,0, 10,  r.)  time  slices.  The  results  in  Fig.  4.5  for  the  Circle  excitation  are 
much  the  same  as  those  for  the  Table  in  Fig.  4.2.  The  "tails,"  however,  shown  in  the  3-D 
image  of  Fig.  4.5  are  rounded  instead  of  cornered  as  in  the  Table  output.  Though  the  Tablt 
output  holds  its  magnitude  for  a  longer  time,  the  maximum  for  the  Circle  output  is  greater 
at  2.21.  (This  value  was  read  from  the  array;  it  is  difficult  to  obtain  quantitative  informa- 
tion from  the  3-d  plots.)  Also  the  drop  off  from  the  maximum  is  steeper  for  the  Circle.  The 
greater  maximum  and  steeper  drop  off  are  due  to  the  equal  distance  of  all  edge  points  from 
the  center.  This  same  geometric  influence  also  accounts  for  the  Circle  holding  the  input 
value  for  a  shorter  duration.  Just  as  the  interference  patterns  added  to  a  maximum  greater 
than  the  input,  the  interference  patterns  also  combine  to  give  a  more  negative  minimum. 
The  negative  undershoot  was  present  for  the  Table;  however,  it  has  a  magnitude  five  times 
greater  for  the  Circle. 

4.3     Other  Input  Excitations 

Having  checked  the  performance  of  the  technique  with  piston  sources,  other  sources  with 
nonuniform  spatial  excitations  were  investigated  The  first  is  a  circularly  truncated  Gaussian 
distribution  function.  Following  the  Gaussian,  a  circularly  truncated  Bessel  profiled  function 
is  examined.  These  two  excitation  function  outputs  are  generated  and  formatted  the  same. 

4.3.1      Gaussian  Distributed  Excitation 

Though  the  Gaussian  has  been  investigated  before  [3]-  [7],  it  has  not  been  studied  as 
a  128x128  array.  The  defining  parameters  for  the  case  studied  are  listed  in  Table  4.4. 
Figure  4.6  shows  the  input  and  resulting  output. 

The  Gaussian  excitation  function  has  been  normalized  by  the  maximum  value  of  the 
computed  Gaussian  (see  the  Gaussian  source  code  titled  CRCGAUS.M  in  Appendix  D). 
This  normalization  is  shown  in  the  input  image  of  Fig.  4.6  by  the  maximum  amplitude  of 
one.  Displayed  as  a  base  x  base  array,  this  input  image  has  a  standard  deviation  of  a  -  5 
and  a  1/e  point  of  10.17  samples  from  the  center  (metric  equivalent  of  r  =  2.54  cm).  The 
diffraction  of  this  input  is  shown  in  Fig.  4.6. 
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(a) 


(b) 


Figure  4.5:  Circle  input  excitation  and  output  for  d  =  31  samples. 
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(o) 


(b) 


Figure  4.6:  Gaussian  distributed  input  and  output  for  a  =  5. 
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NAME 

VALUE 

SUMMARY 

base 

128 

Size  of  square  base  array 

slices 

64 

Number  of  time  slices 

Step 

3 

Number  of  leading  zero-rows 

c 

1500 

Acoustic  velocity  in  media  (m/s) 

z 

0.1 

Distance,  source-to-receiver  (m) 

time.max 

375  x  10-6 

Maximum  time  of  propagation  (s) 

rho.max 

200 

Spatial  radius  of  filters  (length-1) 

d 

51 

Diameter  of  excitation  function  (samples) 

a 

5 

Gaussian  standard  deviation 

Table  4.4:  Defining  parameters  for  Gaussian  excitation  case. 


NAME 

VALUE 

SUMMARY 

base 

128 

Size  of  square  base  array 

slices 

64 

Number  of  time  slices 

Step 

3 

Number  of  leading  zero-rows 

c 

1500 

Acoustic  velocity  in  media  (m/s) 

z 

0.1 

Distance,  source-to-receiver  (m) 

time.max 

375  x  10-6 

Maximum  time  of  propagation  (s) 

rho.max 

200 

Spatial  radius  of  filters  (length-1) 

d 

51 

Diameter  of  excitation  function  (samples) 

a 

0.25 

Bessel  scaling  factor 

Table  4.5:  Defining  parameters  for  Bessel  excitation  case. 


The  3-D  image  shows  a  diffraction  pattern  that  is  well  established  by  time  t  =  time.max, 
forming  two  spreading  "tails."  The  "tails,"  as  well  as  the  rest  of  the  Fig.  4.6.  diffraction 
pattern,  are  smoothly  rounded.  This  rounding  is  the  result  of  the  continuity  of  the  Gaussian 
distribution.  A  discontinuity,  as  in  the  previous  two  excitation  shapes,  results  in  a  char- 
acteristic over  and  undershoot  of  the  maximum  and  minimum  inputs.  Again,  the  results 
for  this  Gaussian  excitation  conformed  to  those  found  in  the  literature  [3,  4,  5,  6,  7].  The 
Bessel  excitation  was  then  run  and  the  results  compared  to  those  of  the  Gaussian. 

4.3.2     Bessel  Excitation 

Results  for  the  Bessel  excitation  produced  by  CRCBES.M  were  generated,  as  previously 
discussed,  for  the  set  of  defining  parameters  listed  in  Table  4.5. 

The  resulting  input  and  output  are  shown  in  Fig.  4.7.     Since  the  output  is  formed  by 
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(a) 


(b) 


Figure  4.7:  Bessel-profile  input  and  output  for  a  =  0.25. 
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taking  successive  p(x,0,10,<)  vectors,  three  peaks  appear  in  the  3-D  image.  These  three 
peaks  correspond  to  the  center  peak  and  the  points  on  the  two  crests  directly  adjacent 
to  the  center.  As  a  result  of  having  three  peaks,  three  "tails"  are  present.  The  "tails," 
however,  are  smooth  because  a  Bessel  function  is  a  continuous  function.  Also  worth  noting 
in  the  3-D  image  is  the  simulated  step  function,  more  visible  here  due  to  the  oscillations  of 
a  Bessel. 

Comparing  the  Gaussian  and  Bessel  outputs,  it  is  seen  that  the  Gaussian's  magnitude 
retention  is  slightly  longer  than  that  of  the  Bessel.  This  is  due  to  the  more  gradual  decrease 
in  magnitude  vice  the  steeper  decrease  required  of  the  Bessel  so  that  it  can  become  negative. 
Still  no  hard  conclusions  can  be  drawn  without  further  analysis. 
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Chapter  5 

Summary 


This  report  presented  a  MATLAB  implementation  of  a  Fourier  approach  to  acoustic  wave 
propagation.  A  mathematical  development  using  linear  systems  that  found  the  acoustic 
potential  from  an  arbitrary  spatial  and  temporal  source  was  reviewed.  In  the  mathematical 
development,  it  was  shown  that  the  Green's  function  solving  the  appropriate  wave  equation 
and  satisfying  the  boundary  conditions  is  the  total  impulse  response  of  the  system.  Through 
double  and  triple  convolutions,  the  acoustic  potential  could  be  found  for  any  source  sepa- 
rable in  time  and  space.  Use  of  the  2-D  spatial  Fourier  transform,  however,  translated  the 
convolution  to  multiplication  in  the  spatial  frequency  domain.  This  made  the  MATLAB 
implementation  easier. 

After  an  overview  of  MATLAB  and  the  graphics  program  AXUM,  a  functional  descrip- 
tion of  the  program  was  furnished.  The  program  modules  AC.FIL  and  AC-PROP  both 
made  use  of  symmetry.  AC-FIL  generated  the  time-varying  filters,  the  most  time  con- 
suming process,  while  AC-PROP  accomplished  the  remaining  computations  making  use  of 
MATLAB's  "fft2"  function.  Details  of  both  modules  as  well  as  the  source  code  have  been 
included  in  the  Appendix  A. 

Several  examples  were  delineated  in  the  body  of  this  thesis.  First,  the  Table  and  the 
Circle  were  presented  as  the  program  verification  excitations;  the  results  conformed  to  those 
found  in  the  literature.  Then  two  newer  excitation  functions,  the  Gaussian  and  the  Bessel, 
were  presented. 

The  underlying  result  was  an  accurate  and  efficient  computer  implementation  of  the 
linear  systems  approach  to  ultrasonic  wave  propagation.  The  efficiency  was  derived  from  the 
modularization  of  the  program  so  that  consecutive  runs  could  be  made  without  recomputing 
the  most  time  consuming  portion,  the  filters.  Also,  the  use  of  MATLAB's  "fft2"  function 
bypassed  tedious  and  time-inefficient  convolution  integrals.  Finally,  both  modules  made 
use  of  symmetry  by  computing  only  one  quadrant  of  data  which  was  then  manipulated  into 
the  remaining  quadrants.  An  advantage  to  using  MATLAB  was  the  ease  of  expansion  that 
could  be  accomplished  with  the  program. 

The  work  of  this  report  concentrated  on  a  source  with  rigidly  baffled  boundary  condi- 
tions and  a  lossless  media.   Cases  that  include  free  space  and  resiliency  baffled  boundary 
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conditions  as  well  as  lossy  media,  linear  and  quadratic  lossy  media,  could  be  incorporated. 
A  few  facts  worth  noting  here  are  that  the  free  space  and  resilient  baffle  boundary  conditions 
can  be  expressed  in  terms  of  the  rigid  baffle  case  [4]  and  that  the  lossy  media  and  lossless 
media  transfer  functions  are  interdependent  [6].  Furthermore,  new  excitation  functions, 
such  as  a  phased  array  or  a  focused  source  [2],  could  also  be  incorporated.  Improvements 
are  needed  in  the  area  of  analysis  such  as  the  Gaussian  versus  Bessel  propagation  compar- 
ison and  extending  the  technique  to  sources  that  are  not  time  and  space  separable  such  as 
new  non-diffracting  waves  [19]. 
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Appendix  A 

Source  Code  for  AC_FIL.M 


The  following  is  the  source  code  used  to  generate  the  time-varying  filters  as  discussed  in 
Chapters  II  and  III.  AC_FIL.M  was  written  in  block  format  with  each  block  headed  by 
a  descriptive  comment  to  explain  the  block's  function.  The  code  includes  many  optional 
instructions  indicated  by  the  leading  %*  symbol.  Deleting  this  symbol  will  enact  that 
line  of  code  on  succeeding  program  runs  which,  in  turn,  varies  the  output.  The  outputs 
necessary  to  a  successful  run  of  AC-PROP.M  are  the  files,  ACbasex(m  +  Step)MAT  (where 
m  is  an  index  number  from  1  to  61).  and  the  file,  AC-FIL.MAT.  For  example,  the  file, 
ACl28x08.MAT,  contains  the  data  for  the  propagation  spatial  filter  in  a  128x128  array  for 
the  fifth  time  slice  (since  Step  =  3).  The  file,  AC-FIL.MAT,  contains  the  parameters  needed 
in  AC.PROP.M. 

Code  is  provided  for  both  the  DOS  version  and  the  Sun  workstation  version  of  MAT- 
LAB4.  There  are  two  primary  differences  in  the  versions.  First,  the  paths  to  disk  storage 
are  different,  to  reflect  the  path  setup  on  the  two  host  machines.  Secondly,  the  commands 
to  obtain  a  hard-copy  version  of  MATLAB4's  graphics  differed.  In  the  DOS  version,  the 
user  prints  copies  from  the  menu  associated  with  the  graphics  window  or  with  a  print  com- 
mand. In  the  SUN  version,  a  hard  copy  can  be  printed  only  by  a  print  command.  Each 
method  is  included  in  the  file  text;  the  user  selects  the  appropriate  command  by  removing 
the  "the  desired  lines. 
ACJPIL.M  SOURCE  CODE 


n 

'/,'/,  This  program  generates  an  Acoustic  Propagation  Transfer 
'/,'/.  Function,  a  time  varying  spatial  filter,  for  use  in 
'/.'/.  AC.PROP.M  to  simulate  acoustic  wave  diffraction. 
VI,  William  H.  Reid        December  1992 

'/.•/,      Modified  for  MATLAB4  and  Sun.  9/93  JPP 

clear;  '/.Clear  MATLAB 
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tml=clock; 


'/.  Start  timer  clock 


base  ■  128; 
NO  =  (base/2)+l; 
slices  =  64; 
Step  =  3; 

c  =  1500; 

2  =  0.1; 

time.max  =  3.75e-4; 

rho.max  =  200; 


Size  of  square  base  array. 

Defines  center  of  square  base  array. 

Number  of  time  slices. 

Number  of  leading  zero  time  slices, 

simulates  the  step  function. 

Velocity  of  the  acoustic  wave,  (m/s) . 

Distance  to  the  observation  plane,  (m) 

Maximum  time  of  propagation, 

time  at  the  final  time  slice,  (sec). 

Spatial  radius  of  the  filter, 

[sqrt(rhox~2  +  rhoy~2] .  (per  length) 


'/,'/,  Initialize  arrays  to  reduce  processing  time, 
shft.filter  =  zeros(base);   temp  =  zeros(NO); 
arg  =  zeros(NO);  rhom  =  zeros(N0,l); 
rho  =  zeros(NO);  time  =  zeros(slices-Step.l) ; 

'/,*/,  Generate  "slices-Step"  time  slices  between  z/c  and  time.max. 
time  ■  linspace(z/c, time.max, slices-Step) ; 

'/,'/,  Choose  directory  to  store  data. 

cd  i:/ac_prop/f ilters  '/,  PC  version 

'/.*  cd  /home2/powers/M_files/ac_prop/f ilters    */,  SUN  version 
'/,'/,  Save  variables  necessary  for  passing  to  AC.PROP.M  in  AC.FIL.MAT. 
save  ac.fil  base  NO  slices  Step  c  z  time.max  rho.max; 

'/,'/,  Generate  N0-1  values  of  "rhom"  from  0  to  rho.max. 
rhom  =  linspace(0,rho_max,N0-l) ; 

'/,'/,  Add  additional  increment  to  rhom  to  compensate  for  off-center 
'/,'/,  orientation  of  the  final  base  x  base  matrix. 
rhom  =  [rhom  (rhom(N0-l)+rhom(2))] ; 

*/,'/,  Create  two  NO  x  NO  arrays  of  rho  values  for  function  evaluation. 
[rhox.rhoy]  ■  meshgrid (rhom, rhom) ; 

VI,   Calculate  "rho",  an  NO  x  NO  matrix  of  radial  distances  for  use  in 
'/,'/,  the  argument  to  the  Bessel  function  within  the  loop. 
rho=  sqrt(rhox."2  +  rhoy.~2); 
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XXXXXXXXXXXX  startloop  xxxxxxxxxx 

XX  Generate  "slices-Step"  filter  arrays,  the  filter  at  each  time  slice. 
for  m  =  1: (slices-Step) 

fprintf  O'/.3.0f '  ,m) ;      '/,  Display  m  value  for  user  progress  report. 

'/,'/,  Create  an  NO  x  NO  array  of  argument  values  for  the  Bessel  function, 
arg  =  rho  *  sqrt(c~2*time(m)~2-z~2  ); 


'/.'/.  Evaluate  the  zero  order  Bessel  for  each  argument  value; 
XX  creates  an  NO  x  NO  array  called  "temp", 
temp  =  2*besselj (O.arg) ; 

VI,   Create  shft.filter  matrix  containing  the  spatial  filter. 
'/,'/,  (Done  by  flipping  "temp"  into  all  quadrants.) 
shft.filter(NO:base,NO:base)  =  tempCl :N0-1 ,1 :N0-1) ; 
shft.filter(NO:base,l:NO-l)  =  f liplr(temp(l :N0-1, 1 :N0-1)) ; 
shft.filter(l:NO-l,l:NO-l)  =  temp(N0:-l :2,N0:-1 :2) ; 
shft.filter(l:NO-l,NO:base)  =  flipud(temp(2:N0,2:N0)) ; 

XX  Save  shft.filter  in  a  file  named  "a(base)x(m+Step) .mat" , 

cd  i:/ac_prop/f ilters  '/,  PC  version 

'/,*  cd  /home2/powers/M_files/ac_prop/f ilters        */,  SUN  version 

if  (m+Step)  <  10, 
•/.'/.  MATLAB  format 

eval(['s3ve  a' ,int2str(base) , 'xO' ,int2str(m+Step) , '  shft.filter  m'  ]  ); 
'/.'/.  ASCII  format 

'/,*  eval(['save  a'  ,int2str(base) ,  'xO'  ,int2str(m+Step) ,  •  .dat  shft.filter. 
'/.*       /ascii']); 

else 
XX  MATLAB  format 

eval(['save  a' ,int2str(base) , 'x' ,int2str(m+Step) , '  shft.filter  m'  ]  ); 
•/.•/.  ASCII  format 

X*  eval(['save  a' ,int2str(base) , 'x' ,int2str(m+Step) , ' .dat  shft.filter.. 
'/.*     /ascii']); 

end 

end 
XXXXXXEND  LOOPXXXXXXXXX 

toe;  X  Stop  timer  clock 

time. elapsed  =  etime(tm2,tml)/60      '/,  Compute  t   display  execution  time 
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Appendix  B 

Examples  of  the  Time— varying 
Propagation  Bessel  Filters 


This  appendix  contains  examples  of  the  time-varying  filters  (Eq.  3.1  on  page  12)  produced 
by  AC_FIL  for  the  parameters  described  in  the  text.  Only  a  few  of  the  61  total  slices  are 
shown.  The  first  slice,  at  t  =  z/c,  is  a  plane  with  amplitude  two;  this  is  because  the  input 
argument  to  the  Bessel  function  is  zero  giving  an  output  value  of  one.  This  uniform  plane 
produces  an  output  for  time  slice  one  that  is  a  scaled  replica  of  the  input.  The  remainder  of 
the  filters  illustrate  the  time  variance  and  how  the  filters  collapse  inward  with  time.  Each 
filter  is  a  128x128  array  representation  in  the  spatial  frequency  domain. 
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(a) 


(b) 


(c) 


Figure  B.l:  Propagation  filter  at  time  slices  1,  2,  and  5. 
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(a 


*£■+ 


(b 


(c) 


Figure  B.2:  Propagation  filter  at  time  slices  10,  20,  and  30. 
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Appendix  C 

Source  Code  for  AC_PROP.M 


The  following  is  the  source  code  for  ACPROP.M  that  was  used  to  produce  various  outputs, 
including  those  discussed  in  Chapter  IV.  AC-PROP  was  written  in  blocks  with  each  block 
headed  with  a  descriptive  comment  to  explain  the  block's  function.  This  also  makes  it 
easy  to  follow  the  computation  from  inputs  to  output.  Inputs  to  AC-PROP  are  imported 
from  the  file  ACFIL.MAT  and  the  files  Abasexm  +  StepM AT  (m  is  an  index  number 
ranging  from  01  to  61).  The  program  solicits  user  input  to  determine  the  input  excitation. 
AC-PROP  then  computes  the  output. 

The  format  of  the  output  can  be  changed  by  the  user.  Before  running  the  program,  the 
user  can  remove  the  optional  comment  markers  indicated  with  %%  to  produce  and/or  view 
the  output  in  the  desired  format.  This  allows  the  data  to  be  saved  and  exported  in  ASCII 
format  for  use  with  other  graphics  programs  (such  as  AXUM,  see  Chapter  II). 

AC-PROP.M  SOURCE  CODE 


,/.,/.,/.,/.,/.  ****  AC.PROP.M  **** 

XX 

'/,'/,  This  script  file  performs  transient-wave  acoustic  propagation 

'/,'/,  simulations.  It  uses  the  time-varying  spatial  filters  called 

'/,'/,  "shft.f ilter"  found  in  the  filters  directory  under  "a(base)x(m)  .mat1 

'/,'/,  to  compute  the  acoustic  propagation  fields.  The  "shft.f ilter" 

'/,'/,  data  are  generated  by  AC.FIL.M. 

XX 

'/.'/,      William  H.  Reid,  December  1992 

•/,'/.      Modified  for  MATLAB4  and  Sun.  9/93  JPP 


clear;  clc;       */,  Clear  MATLAB 
tic;  '/,  Start  timer 


format  compact    */,   Set  compact  format  for  screen  display. 
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'/.'/,  Load  the  parameters  generated  by  last  run  of  AC.FIL.M. 
*/.'/.  PC  and  SUN  directories  have  different  names 

cd  h:\ac_prop\filters  '/,  PC  version 

'/,*  cd  /home2/powers/M_f iles/ac_prop/f ilters  '/,  SUN  version 
load  AC.FIL 

'/,'/,  Display  the  size  of  the  square  base  array. 

dispObase  is  the  width  of  the  array.');  dispC   ');  base 

*/,'/,   Generate  the  INPUT  function  from  user  interface. 

input _func  ■  menu ('Choose  the  shape  of  input  function. ', 'Circle' ,  ... 

'Square' , 'Truncated  Gaussian' , 'Truncated  Bessel'); 
if  isempty(input_func) ;  input.func  =  3;  end  '/,  Default  to  Gaussian  input. 
if  input.func  ==  1,      '/,  Cicle  input 
name='c' ; 

dispC'You  have  chosen  a  truncated  circle  input.') 
d  =  input ('Please  enter  an  ODD  diameter,  [51],  d  =   '); 

if  isempty(d);  d  =  51;  end  '/,  Default  diameter:  51  samples 
shft.input  =  crcle(d.base) ;    */.  Create  Circle  input 
'/,  Name  a  file  to  hold  the  input  function. 
InFile.Name  =  [name, 'i' ,int2str(base) , 'x' ,int2str(d)] ; 
'/,  Name  a  file  to  hold  the  output. 

OutFile.Name  =  [name, 'o' ,int2str(base) , 'x' ,int2str(d)] ; 
elseif  input.func  ==  2,   '/,  Square  input 
name='s' ; 

dispCYou  have  chosen  a  truncated  square  input.') 
v  =  input ('Please  enter  an  ODD  width,  [11],  w  =   '); 

if  isempty(w);  w  *  11;  end  '/,  Default  width:  11  samples 
shft.input  =  table (w, base) ;    '/,  Create  input 

d  =  w; 
'/,  Name  a  file  to  hold  the  input  function. 
InFile.Name  =  [name, 'i' ,int2str(base) , 'x' ,int2str(d)] ; 
'/,  Name  a  file  to  hold  the  output. 

OutFile_Name  =  [name, 'o' ,int2str(base) , 'x' ,int2str(d)] ; 
elseif  input.func  ■■  3,    '/,  Gaussian  input 
name='g' ; 

dispCYou  have  chosen  a  truncated  Gaussian  input.') 
sigma  =  input('Please  enter  the  standard  deviation,  [10],  sigma  =  '); 

if  isempty (sigma) ;  sigma  =  10;   end  '/.  Default  sigma:  10  samples 
d  ■  input ('Please  enter  an  ODD  diameter,  [51],  d  =   '); 

if  isempty (d);  d  =  51;  end  '/,  Default  diameter:  51  samples 
shft.input  =  crcgaus(sigma,d,base) ;   */,  Calculate  input 
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'/,  Name  a  file  to  hold  the  input  function. 

InFile.Name  =  [name, 'i' ,int2str(base) , *x' ,int2str(d)] ; 

'/,  Name  a  file  to  hold  the  output. 

OutFile.Name  =  [name, 'o' ,int2str(base) , 'x' ,int2str(d)] ; 
elseif  input.func  ==  4,     '/,  Bessel  input 

name='b' 

disp('You  have  chosen  a  truncated  Bessel  input.') 

a  =  input('Please  enter  a  width  scaling  factor, [0.3125] ,  a  =   '); 
if  isempty(a);  a  =  0.3125;  end  '/.  Default:  0.3125 

d  =  input('Please  enter  the  ODD  diameter,  [51],  d  =   '); 
if  isempty(d);  d  =  51;  end  '/.  Default:  51  samples 

disp( 'Please  wait.  This  calculation  takes  a  while ') 

shft.input  =  crcbess(a,d,base) ;  */.  Calculate  input 

q  =  a  *  le4; 

'/,  Name  a  file  to  hold  the  input  function. 

InFile.Name  =  [name, 'i' ,int2str(base) , 'x' ,int2str(q)] ; 

'/,  Name  a  file  to  hold  the  output. 

OutFile.Name  =  [name, 'o' ,int2str(base) , 'x' ,int2str(q)]  ; 
lse 

dispC    >); 

disp( 'Incorrect  Excitation  Function  Selection!'); 

error ('Restart  AC. PROP... to  try  again.'); 
end 

'/.'/,  Save  the  input  function  placed  in  "InFile.Name"  in  ascii  format 
'/,'/,  for  use  with  other  graphics  software. 

eval(['save  ' .InFile.Name, ' .dat  shft.input  /ascii']); 

clear  InFile.Name;       '/,  Remove  "InFile.Name"  from  memory. 

'/,'/,  Compute  the  sample  spacing  in  x  and  y  (delta.x)  and  the  spacing  in 
VI,         time  (delta.t) . 

delta.x  ■  (base-2)/(2*rho_max) ;  delta.t  =  time.max/ (slices-Step) ; 
'/,'/,  Create  X  and  T  vectors  (1  x  base). 

X  =  linspace(-(NO-l)*delta_x, ((base-NO)-l)*delta_x,base) ; 

T  =  linspace(-3*delta_t, time.max, slices) ; 

*/,'/,  Display  the  input  function. 

'/.*  mesh(X,X,abs(shft_input));title('|SHFT.INPUT|'); 

*/.*  axis(  [  X(l)  X(base)  X(l)  X(base)  ... 

'/,*  min(    [  0  min(min(abs (shft.input)))   ]    )  max(max(abs (shft.input)))   ]) 

*/.*  xlabeK'x   (cm)');   ylabeK'y   (cm)');   zlabeK  'amplitude' ) ;   view(52.5,30) ; 
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'/,'/,  SUN:  Save  input  figure  as  eps  file. 

'/,*  cd  /home2/powers/M_f iles/ac.prop/data 

'/,'/,  PC:  save  figure  from  figure  window. 

'/,*  cd  h:/ac.prop/ 

'/,'/.  Save  input  plot  as  EPSF  file  with  bitmap  preview  image. 

'/.*  dispC  ');  dispOSaving  input  plot  as  EPS  file ');  dispC  '); 

*/,*  eval(['print  -deps  -epsi  '  ,name,  '_in_b'  ,int2str(base) ,  ... 
'/.*  '.d',int2str(d)]); 

*/,*  dispC  ');  dispCSaving  input  plot  as  EPS  file ');  dispC  '); 

'/,*  eval(['print  -deps  -epsi  '  .name, '_in_b' ,int2str(base) ,  ... 
'/.*  '_d',int2str(d)]); 

'/,'/,  Shift  "shft.input"  from  centered  geometry  to  corner  geometry  and  take 
VI,         2-D  FFT  to  produce  F_ INPUT. 

F.input  =  fft2(  fftshift(shft.input)  ); 

clear  shft.input;  '/,  Free  RAM. 

'/,'/,  Shift  F.input  in  preparation  of  multiplication  with  PROP.m. 
shft.F.input  =  fftshift (F.input) ; 
clear  F.input;   '/.  Free  RAM. 

'/.'/.  Element  by  element  arrray  multiplication  of  the  transfer  function 
'/.'/.  filter  in  "PROP.m"  and  "Fshft. input .  " 
dispC  Performing  array  multiplication....'); 


'/,*  cd  /home2/powers/M_f iles/ac.prop/data  */,  SUN  version 
cd  h:/ac_prop/data  '/,  PC  version 

VI,   Save  "Step"  arrays  full  of  zeros  for  first  array  values . 

vi,vi:i:i:i:i:i:i,vi:i,  start  first  loop  vaivi,vi:i:i:i:i:i:a 

for  m=l :Step 

shft.output  =  zeros (base); 
VI,   Save  the  (base)x(base)  array  "shft.output"  in  an  ASCII  file  whose 
•/,'/,  name  is  "(input_func)_OUT_(m)  .dat"  . 

eval(['save  ' ,int2str(input_func) , '.OUT.O' ,int2str(m) ,  ... 
'.dat  shft.output  /ascii']); 

clear  shft.output;       '/,  Get  ready  for  next  pass, 
end 

vi.vi:i,vi:i:i:i:i:i:i.  end  first  loop  vi,vi,vi,vi:i:i:i:i:i:i:i, 
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xmmmxxx  start  second  loop  vm:m:m:i:i:i:ii 

for  m  =  1 : (slices-Step) 

fprintf  (,,/,2.0f , '  ,m)      '/,  Display  "m"  value  for  progress  report. 

VI,   Give  the  variable  "filenamel"  the  name  of  the  file  containing 
XX  "shft.filter"  and  then  load  that  file, 
if  m  <  10-Step, 

filenamel  =  ['a' ,int2str(base) , 'xO' ,int2str(m+Step)] ; 
else 

filenamel  =  ['a' ,int2str(base) , 'x' ,int2str(m+Step)] ; 
end 
*/,*  cd  /home2/powers/M_files/ac_prop/f ilters  '/,  SUN  version 
cd  h:/ac_prop/f ilters  '/,  PC  version 

eval(['load  • .filenamel]  ); 

*/,'/,  Multiply  the  filter  by  the  shifted  transform  of  the  input. 
shft.F.output  =  (shft.filter  .*  (shft.F.input)) ; 
clear  shft.filter;  X  Free  RAM 

'/.'/,  Shift  "shft.F.output"  from  centerd  geometry  to  corner  geometry 
XX    ("F.output")  and  take  inverse  transform  to  produce  "output." 

output  =  ifft2(  fftshift (shft.F.output)  )  ; 
'/.'/,  Shift  "output"  to  center  geometry  ("shft.output") . 

shft.output  =  fftshift(output) ; 
'/.'/,  "shft.output"  is  the  centered  geometry  version  of  the  diffracted 
'/.'/,      wave  at  time  slice  "m" . 

clear  shft.F.output  output  '/,  Free  RAM. 

XX  Place  the  transposed  "NO"  (center)  row  of  "shft.output"  in  the 
XX  (m+Step)-th  column  of  "output .plot . " 

'/,'/,  This  creates  a  (base)  x  (m+Step)  array  whose  columns  show  a  slice 
*/.'/,  of  the  diffracted  wave. 

output .plot (l:base,m+ Step)  =  (shft.output (NO, 1 :base)) ' ; 

XX  Save  the  (base)x(base)  array  "shft.output"  in  an  ASCII  file  whose  name  is 
•/•/.  "  (input _func) .OUT. (m+Step)  .dat". 

XX  A  Gaussian,  for  example,  would  be  G_0UT_12.dat  on  the  9-th  loop 
VI,   interation  (Step  =  3) .  This  is  optional  for  graphics  use  by 
'/.'/,  other  software. 

X*  cd  /home2/powers/M_f iles/ac.prop/data  '/,  SUN  version 
cd  h:/ac_prop/data  '/,  PC  version 

if  m  <  10-Step, 
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eval(['save  ' ,int2str(input_func) , '.OUT.O' ,int2str(m+Step) ,  ... 
'.dat  shft_output  /ascii']); 
else 

eval(['save  ' ,int2str(input_func) , '.OUT.' ,int2str(m+Step) ,  ... 
'.dat  shft.output  /ascii']); 
end 
clear  shft.output;       X  Get  ready  for  next  pass, 
end 

xxxxxxxxxxxxxxxxx  End  iooP  xxxxxxxxxxxxxxx 

*/,'/,  Save  the  contents  of  "output.plot"  in  a  MATLAB  file 

XX  and  as  an  ascii  file  (optional) . 

'/,*  cd  /home2/powers/M_f  iles/ac.prop/data  X  SUN  version 
cd  h:/ac_prop/data  '/,  PC  version 

eval(['save  '.OutFile.Name,'  output.plot']) 
eval(['save  ' .OutFile.Name, ' .dat  output.plot  /ascii']); 

'/,'/,  Display  output  view  #1. 

XX  figure;  mesh(T, X, abs (output.plot) ) ;  titleC I0UTPUTI ') ; 

XX  axis(  [  T(l)  T(slices)  X(l)  X(base)   ... 

XX  min(  [  0  min(min(abs (output.plot)))  ]  )  ... 

XX  max (max (abs (output.plot)))  ]) 

XX  xlabeK'time   (s)');   ylabeK'y   (cm)');   zlabeK 'abs  (output) ') ; 

XX  view (52. 5, 30); 

XX  Save  output  figure  as  eps  file. 

X*  cd  /home2/powers/M_f iles/ac.prop/data  X  SUN  version 
cd  h:/ac_prop/data  X  PC  version 

X*  dispC  ');  disp( 'Saving  output  figure  as  EPS  file ');  dispC  '); 

X*  eval(['print  -deps  -epsi  ' .name, '.outl.b' ,int2str (base) ,  ... 
X*  '_d\  int2str(d)]); 

XX  Create  a  different  view  of  output 

X*   figure;  mesh(T,X,abs(output.plot)) ;  titleC I0UTPUTI ') ; 

X*   axis(  [  T(l)  T(slices)  X(l)  X(base)   ... 

X*       min(  [  0  min(min(abs (output.plot)))  ]  )  ... 

X*       max (max (abs (output.plot)))  ]) 

X*  xlabeK'time  (s)');  ylabeK'y  (cm)');  zlabel  Cabs  (output) ') ; 

XX  Save  second  output  figure  as  eps  file. 

X*  dispC  *);  dispC  Saving  figure  as  EPSF  file ');  dispC  '); 

X*  evaKC'print  -deps  -epsi  '  .name,  '_out2_b'  ,int2str(base) ,  '_d' ,  ... 
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'/.*      int2str(d)]); 

elapsed.time  =  toc/60  '/,  Stop  timer;  display  elapsed  time 
disp( 'minutes') 
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Appendix  D 

Source  Code  for  Input  excitations 


The  following  source  code  is  for  the  input  excitation  choices  given  to  the  user  in  AC-PROP. 
Each  is  written  as  a  MATLAB  function  and  can  be  used  independently  of  AC.FIL  or 
AC-PROP.  The  input  excitations  are  the  uniform  square  (TABLE(w,base)),  the  uniform 
circle  (CRCLE(d,base)),  the  circularly  truncated  Gaussian  (CRCGAUS(sigma,d,base)),  and 
the  circularly  truncated  Bessel  (CRCBES(a,d,base))  where  w  is  the  width  of  the  square,  d 
is  the  diameter  of  the  circle,  a  is  the  standard  deviation  of  the  Gaussian  distribution,  a  is 
a  scaling  factor,  and  base  is  the  size  of  the  base  array. 

TABLE.M  SOURCE  CODE 


function  Y  =  table(w,base) 

V.     TABLE.M:  Y  =  table(w.base) 

'/.Program  for  generating  a  uniform  square  excitation  function. 

'/,   December  1992     William  H.  Reid 

'/,   Based  on  TABLE.M  by  JG  Upton 

y. 

X   w  is  the  WIDTH  of  the  table.   (ODD  integer) 

'/,   base  is  the  WIDTH  of  the  square  base.   (EVEN  integer) 

'/.   Example:  z  =  table(33,64) ; 

'/,  Check  that  w  is  an  odd  integer, 
if  rem(w,2)  <  0.1; 

error('The  width  of  the  table  must  be  an  ODD  integer.'); 
else; 
end; 

'/,  Check  that  base  is  an  even  integer. 
if  rem(base,2)  "■  0.0; 
errorCThe  width  of  the  square  base  must  be  an  EVEN  integer.'); 
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else; 
end; 

NO  =   (base/2)+l;  '/,  NO  is  the  base  array's   center. 

wO  =  ceil(w/2);  '/,  wO  is  the  mid--point  of  the  table 

Y  =  zeros(base);  '/,  Initialize  matrices  to  reduce 

temp  =  zeros(NO-l);  '/,  processing  time. 

temp(l:wO,l:uO)   =  ones(wO);  '/,  Set  amplitude  to  one. 

'/,  Generate  the  entire  base  x  base  input  function  array. 
Y(baseO:base,baseO:base)   =  temp; 
Y(2:base0,base0:base)   =  rot90(temp); 
Y(2:N0,2:N0)   =  rot90(temp,2) ; 
Y(N0:base,2:N0)   =  rot90(temp,3) ; 

'/,  To  test   input  distribution:     mesh(Y) 


CRCLE.M  SOURCE  CODE 


function  Y  =  crcle(d.base) 

CRCLE.M:   Y  =  crcle(d.base) 
Program  for  generating  uniform  circular  excitation  functions 
December  1992     William  H.  Reid 
Based  on  CRCLE.M  by  JG  Upton 

d  is  the  DIAMETER  of  the  circle.   (ODD  integer) 
base  is  the  WIDTH  of  the  square  base.  (EVEN  integer) 
Example:  z  =  crcle(33,64) ; 

'/.  Check  that  d  is  an  odd  integer, 
if  rem(d,2)  <  0.1; 

error ('The  diameter  of  crcle  function  must  be  an  ODD  integer.'); 
else; 
end; 

'/,  Check  that  base  is  an  even  integer. 
if  rem(base,2)  "=  0.0; 
error('The  width  of  the  square  base  must  be  an  EVEN  integer.'); 
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else; 
end; 

NO  =  (base/2)  +  l;  */,  NO  is  the  base  array's  center. 

r  =  d/2;  '/,  r  is  the  circle's  radius. 

'/,  Initialize  matrices  to  reduce  processing  time. 
Y  =  zeros(base) ; 
temp  =  zeros(NO-l); 

'/,  Set  amplitude  to  one  inside  the  circle's  radius. 
for  m  =  1 :r+l ; 

for  n  =  l:r+l; 

if  sqrt((m-l)~2  +  (n-l)"2)  <=  r; 
temp(m.n)  =  1; 
end; 
end; 
end; 

'/,  Generate  the   entire  base  x  base   input  function. 
Y(baseO:base ,baseO :base)    =  temp; 
Y(2:base0,base0:base)    =  f lipud(temp) ; 
Y(2:N0,2:N0)    =  rot90(temp,2) ; 
Y(N0:base,2:N0)    =  fliplr(temp) ; 

7,  To  test   input  function  distribution:     mesh(Y) 

CRCGAUS.M  SOURCE  CODE 

unction  Y  =   crcgaus(sigma,d,base) 

CRCGAUS.M:      Y  =  crcgaus(sigma,d,base) 
Program  for  generating  circular  Gaussian  excitation  functions. 
December   1992  William  H.   Reid 

Based  on  CRCGAUS.M  by  JG  Upton 

sigma  is  the  STANDARD  DEVIATION  of  the  Gaussian  function. 
d  is  the     DIAMETER  of   circle.      (ODD  integer) 
base  is  the  WIDTH  of  the  square  base.    (EVEN   integer) 
Example:   z  ■  crcgaus(12,33,64) ; 

mu=0;  y.mu  is  the  mean  of  the  Gaussian  function 
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'/.  Check  that  d  is  an  odd  integer, 
if  rem(d,2)  <  0.1; 

error('The  diameter  of  circle  function  must  be  an  ODD  integer.'); 
else; 
end; 

'/,  Check  that  base  is  an  even  integer, 
if  rem(base,2)  "=  0.0; 

error('The  width  of  the  square  base  must  be  an  EVEN  integer.'); 
else; 
end; 

NO  =  (base/2)  +  l;  */.  NO  is  center  of  the  array. 

r  =  d/2;  '/,  r  is  the  radius  of  the  truncating  circle 


*/,  Initialize  the  matrices  to  reduce  processing  time. 

Y  =  zeros (base) ; 
temp  =  zeros(NO-l); 

*/,  Compute  the  amplitude  for  the  Gaussian  distributed  circle, 
for  m  =  l:(d+l)/2; 

for  n  =  l:(d+l)/2; 

x  =  sqrt((m-l)*2+(n-l)*2); 

if  x  <=  r; 

temp(m,n)  =  (l/(sqrt(2*pi)*sigma))*exp(-((x-mu)~2)/. . . 

(2*(sigma~2))); 
end; 
end; 
end; 

'/,  Generate  the  entire  base  x  base  input  array. 
Y(baseO:base,baseO:base)  =  temp; 
Y(2:base0,base0:base)  =  flipud(temp) ; 
Y(2:N0,2:N0)  =  rot90(temp,2) ; 
Y(N0:base,2:N0)  =  fliplr(temp) ; 

Y  =  Y  ./  (max(max(Y))) ;  */,  Normalize  the  Gaussian  distribution  to  one 
'/.  To  test  and  view  the  input  function:  mesh(Y) 
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CRCBESS.M  SOURCE  CODE 


function  Y  =  crcbess(a,d,base) 

CRCBESS.M:   Y  =  crcbess(a,d,base) 
Program  for  generating  circular  Bessel  excitation  functions. 
December  1992     William  H.  Reid 
Based  on  CRCBESS.M  by  JG  Upton 

a  is  the  WIDTH  SCALING  FACTOR. 
d  is  the  DIAMETER  of  the  circle.   (ODD  integer) 
base  is  the  WIDTH  of  the  square  base.  (EVEN  integer) 
Example:  z  =  crcbessd  ,33,64) ; 

'/,  Check  that  d  is  an  odd  integer. 
if  rem(d,2)  <  0.1; 

error ('The  diameter  of  the  circle  must  be  an  ODD  integer'); 
else; 
end; 

'/,  Check  that  base  is  an  even  integer. 
if  rem(base,2)  "=  0.0; 

errorCThe  width  of  the  square  base  must  be  an  EVEN  integer'); 
else; 
end; 

NO  =  (base/2)  +  l;  '/,  NO  is  the  center  of  the  array. 

r  =  d/2;  '/,  r  is  the  radius  of  the  circle. 

Y  =  zeros(base);  '/,  Initialize  the  arrays  to  reduce 

temp  =  zeros (N0-1);  '/,  processing  time. 

'/.  Compute  the  Bessel  distributed  amplitude  within  the  circle, 
for  m  =  l:r+l; 

for  n  =  l:r+l; 

x  =  sqrt((m-l)"2  +  (n-l)~2); 
if  x  <=  r; 

temp(m,n)=besseln(0,a*x) ; 
end; 
end; 
end; 

'/,  Generate  the  entire  base  x  base  input  array. 
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Y(NO:base,NO:base)  =  temp; 
Y(2:N0,N0:base)  =  flipud(temp) ; 
Y(2:N0,2:N0)  =  rot90(temp,2) ; 
Y(N0:base,2:N0)  =  fliplr(temp) ; 

'/,  To  test  and  view  the  input  function:  mesh(Y) 
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Appendix  E 

Examples  of  Dicer 
Representations  of  Output  Data 


This  appendix  contains  more  representations  of  the  data  from  the  Spyglass  Dicer  program. 
The  representations  are  discussed  in  the  text  on  page  4.3. 

Figures  E.l  and  E.2  show  representations  of  data  from  a  Table  excitation  that  is  23 
samples  wide.  Figures  E.3  and  E.4  show  representations  of  data  from  a  Circle  excitation 
that  is  23  samples  wide.  Figures  E.5  and  E.6  show  representations  of  data  from  a  Circle 
excitation  that  is  31  samples  wide. 
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Figure  E.l:  Table  output  for  d  =  23  samples.  Dicer  representation. 
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Figure  E.2:  Table  output  for  d  =  23  samples.  Alternative  Dicer  representation. 
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Figure  E.3:  Circle  output  for  d  -  23  samples.  Dicer  representation. 
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Figure  E.4:  Circle  output  for  d  =  23  samples.  Alternative  Dicer  representation. 
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Figure  E.5:  Circle  output  for  d  =  31  samples.  Dicer  representation. 
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Figure  E.6:  Circle  output  for  d  =  31  samples.  Alternative  Dicer  representation. 
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